################################################################################################%
# This script takes the bootstrapped SV data under the rules specified in the SV Welfare  
# Simulations Script and then produces inequality measures and graphs.
# 
# All graphs are exported to the SV & QV Recalibration on Pop Imm/SV - Inequality/ .
# 
# Author: Luis S.
# Last Modified: 6.19.18
################################################################################################%

library(dplyr)
library(ineq)
library(ggplot2)
library(stargazer)

setwd("C:/Users/Luis/Dropbox/Research/SV vs QV California Experiment/Experiment Analysis & Results")


SV_IndivUtility_NSV_A <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_A.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_NSV_B <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_B.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_NSV_C <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_C.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_NSV_D <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_D.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_NSV_E <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_E.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV1_A <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_A.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV2_A <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_A.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV1_B <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_B.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV2_B <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_B.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV1_C <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_C.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV2_C <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_C.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV1_D <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_D.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV2_D <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_D.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV1_E <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_E.csv", header = T, stringsAsFactors = F)
SV_IndivUtility_SV2_E <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_E.csv", header = T, stringsAsFactors = F)

SV_RuleAResults <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleAResults.csv", header = T, stringsAsFactors = F)
SV_RuleBResults <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleBResults.csv", header = T, stringsAsFactors = F)
SV_RuleCResults <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleCResults.csv", header = T, stringsAsFactors = F)
SV_RuleDResults <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleDResults.csv", header = T, stringsAsFactors = F)
SV_RuleEResults <- read.csv("SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleEResults.csv", header = T, stringsAsFactors = F)

SV_RuleAResults <- mutate(SV_RuleAResults, Rule = "A")
SV_RuleBResults <- mutate(SV_RuleBResults, Rule = "B")
SV_RuleCResults <- mutate(SV_RuleCResults, Rule = "C")
SV_RuleDResults <- mutate(SV_RuleDResults, Rule = "D")
SV_RuleEResults <- mutate(SV_RuleEResults, Rule = "E")

SV_AllRules <- bind_rows(SV_RuleAResults, SV_RuleBResults, SV_RuleCResults, SV_RuleDResults, SV_RuleEResults)

####################################
###### Determine Minority Win ######
###################################%

SV_AllRules <- SV_AllRules %>% 
  mutate(SV1_MinorityWon = ifelse(SV1_Winner_BilingualEduc == "Minority" | SV1_Winner_Immigration == "Minority" | 
                                    SV1_Winner_PublicBonds == "Minority" | SV1_Winner_TeacherTenure == "Minority",T, F),
         SV2_MinorityWon = ifelse(SV2_Winner_BilingualEduc == "Minority" | SV2_Winner_Immigration == "Minority" | 
                                    SV2_Winner_PublicBonds == "Minority" | SV2_Winner_TeacherTenure == "Minority",T, F))


SV_RuleAResults <- filter(SV_AllRules, Rule == "A")
SV_RuleBResults <- filter(SV_AllRules, Rule == "B")
SV_RuleCResults <- filter(SV_AllRules, Rule == "C")
SV_RuleDResults <- filter(SV_AllRules, Rule == "D")
SV_RuleEResults <- filter(SV_AllRules, Rule == "E")

bootstrapIterations <- nrow(SV_RuleAResults)


###################################################
###### Inequality Measure - Gini Coefficient ######
##################################################%

### Testing
# sampleData <- SV_AllRules
# rule <- "B"
# vote <- "SV1"
# util_NSV <- SV_IndivUtility_NSV_B
# util_SV <- SV_IndivUtility_SV1_B
 

getGiniCoeff <- function(sampleData, rule, vote, util_NSV, util_SV)
{
  # Note we use PercentGiniChange = (Gini_UtilSV - Gini_UtilNSV)/Gini_UtilNSV
  tempRule <- filter(sampleData, Rule == rule)
  
  # get list of simulations in which minority won
  roundsWithMinVic <- tempRule[tempRule[[paste(vote,"_MinorityWon", sep = "")]] == TRUE,]
  roundsWithMinVic <- roundsWithMinVic %>% select(Round) %>% mutate(MinVic = T)
  
  # get utilities under SV and Simple Maj when minority won with SV
  tempRule_UtilNSV <- inner_join(select(tempRule, Round), util_NSV, by = "Round") %>% select(-matches("Round"))
  tempRule_UtilSV <- inner_join(select(tempRule, Round), util_SV, by = "Round") %>% select(-matches("Round"))
  
  # create dataframe to store gini coeff
  tempRule_GiniCoeffs <- data.frame(Round = tempRule$Round, Rule = rule, Gini_UtilNSV = numeric(nrow(tempRule)), Gini_UtilSV = numeric(nrow(tempRule)))
  
  # iterate over each simulation in which minority won to get gini coeff
  for(row in 1:nrow(tempRule))
  {
    print(row)
    tempRule_GiniCoeffs$Gini_UtilNSV[row] <- Gini(tempRule_UtilNSV[row,], corr = FALSE)
    tempRule_GiniCoeffs$Gini_UtilSV[row] <- Gini(tempRule_UtilSV[row,], corr = FALSE)
  }
  
  # Join with RoundsWithMinVic 
  tempRule_GiniCoeffs <- full_join(tempRule_GiniCoeffs, roundsWithMinVic, by = "Round")
  tempRule_GiniCoeffs[is.na(tempRule_GiniCoeffs)] <- F 
  
  # Calculate percentage Gini change between NSV and SV1 & and determine again whether Gini decreased or not
  tempRule_GiniCoeffs <- mutate(tempRule_GiniCoeffs, Gini_UtilNSV = round(Gini_UtilNSV, digits = 4), Gini_UtilSV = round(Gini_UtilSV, digits = 4)) %>%
    mutate(PercentGiniChange = (Gini_UtilSV - Gini_UtilNSV)/Gini_UtilNSV, PercentGiniChange = round(PercentGiniChange, digits = 4),
           ChangeDirection = ifelse(Gini_UtilSV == Gini_UtilNSV, "No Change", ifelse(Gini_UtilSV < Gini_UtilNSV, "Decrease", "Increase")),
           GiniDecrease = ifelse(ChangeDirection == "Decrease", TRUE, FALSE))
  
  
  # Save Gini Raw Data 
  write.csv(tempRule_GiniCoeffs, paste("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - ",vote,"_Rule",rule,"_GiniCoeffs.csv", sep = ""), row.names = F)
  
  return(tempRule_GiniCoeffs)
}


#--- SV1: Calculating Gini Coefficient, percent gini change, & freq of gini decline 

# Calculate Gini coeffs for all sims and all rules

SV1_RuleA_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "A", vote = "SV1", util_NSV = SV_IndivUtility_NSV_A, util_SV = SV_IndivUtility_SV1_A) 
SV1_RuleB_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "B", vote = "SV1", util_NSV = SV_IndivUtility_NSV_B, util_SV = SV_IndivUtility_SV1_B)
SV1_RuleC_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "C", vote = "SV1", util_NSV = SV_IndivUtility_NSV_C, util_SV = SV_IndivUtility_SV1_C)
SV1_RuleD_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "D", vote = "SV1", util_NSV = SV_IndivUtility_NSV_D, util_SV = SV_IndivUtility_SV1_D)
SV1_RuleE_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "E", vote = "SV1", util_NSV = SV_IndivUtility_NSV_E, util_SV = SV_IndivUtility_SV1_E)


#--- Modified ---#

# SV1_RuleA_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV1_RuleA_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# SV1_RuleB_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV1_RuleB_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# SV1_RuleC_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV1_RuleC_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# SV1_RuleD_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV1_RuleD_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# SV1_RuleE_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV1_RuleE_GiniCoeffs.csv", stringsAsFactors = F, header = T)

meanGiniSummary <- data.frame()

getNthEmpPercentile <- function(nums, NthPercent)
{
  orderedNums <- sort(nums)
  totalNums <- length(orderedNums)
  
  closestPercentile <- round((totalNums/100)*NthPercent)
  
  return(orderedNums[closestPercentile])
}

# Gini of Average Realized Utilities (over all 10k sims)

SV1_AllRules_GiniCoeffs <- bind_rows(SV1_RuleA_GiniCoeffs, SV1_RuleB_GiniCoeffs, SV1_RuleC_GiniCoeffs, SV1_RuleD_GiniCoeffs, SV1_RuleE_GiniCoeffs)

meanGiniSummary <- bind_rows(meanGiniSummary, SV1_AllRules_GiniCoeffs %>% group_by(Rule) %>% 
                               summarize(MeanGini_Maj = mean(Gini_UtilNSV), lower95CI_Maj = getNthEmpPercentile(Gini_UtilNSV, 2.5), upper95CI_Maj = getNthEmpPercentile(Gini_UtilNSV, 97.5),
                                         MeanGini_SV = mean(Gini_UtilSV), lower95CI_SV = getNthEmpPercentile(Gini_UtilSV, 2.5), upper95CI_SV = getNthEmpPercentile(Gini_UtilSV, 97.5)) %>%
                               mutate(Sample = "SV1", AveragedOver = "All Sims"))


# Filter data so we only make plots conditional on a minority victory

SV1_RuleA_GiniCoeffs <- SV1_RuleA_GiniCoeffs %>% filter(MinVic == T) 
SV1_RuleB_GiniCoeffs <- SV1_RuleB_GiniCoeffs %>% filter(MinVic == T) 
SV1_RuleC_GiniCoeffs <- SV1_RuleC_GiniCoeffs %>% filter(MinVic == T) 
SV1_RuleD_GiniCoeffs <- SV1_RuleD_GiniCoeffs %>% filter(MinVic == T) 
SV1_RuleE_GiniCoeffs <- SV1_RuleE_GiniCoeffs %>% filter(MinVic == T) 


# Gini of Average Realized Utilities (over sims with min vic)

SV1_AllRules_GiniCoeffs <- bind_rows(SV1_RuleA_GiniCoeffs, SV1_RuleB_GiniCoeffs, SV1_RuleC_GiniCoeffs, SV1_RuleD_GiniCoeffs, SV1_RuleE_GiniCoeffs)

meanGiniSummary <- bind_rows(meanGiniSummary, SV1_AllRules_GiniCoeffs %>% group_by(Rule) %>% 
                               summarize(MeanGini_Maj = mean(Gini_UtilNSV), lower95CI_Maj = getNthEmpPercentile(Gini_UtilNSV, 2.5), upper95CI_Maj = getNthEmpPercentile(Gini_UtilNSV, 97.5),
                                         MeanGini_SV = mean(Gini_UtilSV), lower95CI_SV = getNthEmpPercentile(Gini_UtilSV, 2.5), upper95CI_SV = getNthEmpPercentile(Gini_UtilSV, 97.5)) %>%
                               mutate(Sample = "SV1", AveragedOver = "Miv Vic"))


#--- Modified (end) ---#


# Filter data so we only make plots conditional on a minority victory

SV1_RuleA_GiniCoeffs <- SV1_RuleA_GiniCoeffs %>% filter(MinVic == T) 
SV1_RuleB_GiniCoeffs <- SV1_RuleB_GiniCoeffs %>% filter(MinVic == T) 
SV1_RuleC_GiniCoeffs <- SV1_RuleC_GiniCoeffs %>% filter(MinVic == T) 
SV1_RuleD_GiniCoeffs <- SV1_RuleD_GiniCoeffs %>% filter(MinVic == T) 
SV1_RuleE_GiniCoeffs <- SV1_RuleE_GiniCoeffs %>% filter(MinVic == T) 

# SV1 - Summary of Gini Changes

SV1_AllRules_GiniCoeffs <- bind_rows(SV1_RuleA_GiniCoeffs, SV1_RuleB_GiniCoeffs, SV1_RuleC_GiniCoeffs, SV1_RuleD_GiniCoeffs, SV1_RuleE_GiniCoeffs)

rulesSummary <- SV1_AllRules_GiniCoeffs %>% group_by(Rule, ChangeDirection) %>% summarize(Count = n()) %>% mutate(Prop = round(Count/sum(Count), digits = 4))

stargazer(rulesSummary, 
          title = "Summary of Gini Changes by Rule (SV1)", type = "text", digits = 3, summary = FALSE, rownames = FALSE, 
          out = "./SV & QV Recalibration on Pop Imm/SV - Inequality/Table - Summary of Gini Changes by Rule (SV1).txt")

# SV1 - Dot Plot of Frequency of Gini Declines per Rule in (Fig 8 - top) 
FreqGiniDeclineA <- group_by(SV1_RuleA_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineB <- group_by(SV1_RuleB_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineC <- group_by(SV1_RuleC_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineD <- group_by(SV1_RuleD_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineE <- group_by(SV1_RuleE_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)

SV1_FreqGiniDeclines <- data.frame(Rule = c("A", "B", "C", "D", "E"), 
                                   Freq = c(FreqGiniDeclineA$Freq[1], FreqGiniDeclineB$Freq[1], FreqGiniDeclineC$Freq[1], 
                                            FreqGiniDeclineD$Freq[1], FreqGiniDeclineE$Freq[1]))
SV1_FreqGiniDeclines <- mutate(SV1_FreqGiniDeclines, Freq = round(Freq, digits = 3)*100)

inequalityPlot <- ggplot(data = SV1_FreqGiniDeclines, aes(x=Rule, y = Freq)) + geom_point() + 
  geom_text(aes(label=Freq), hjust = .5, vjust = 2) + ggtitle("SV1 - Percent of Gini Declines by Voting Rule") + 
  scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Freq. of Gini Declines (%)") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12))
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/SV - Inequality/SV1 - Percent of Gini Declines by Voting Rule.png", height = 5.5, width = 5.5)


# SV1 - Bar Plot of Mean Percent Gini Change per Rule in (Fig 8 - bottom) 
SV1_MeanGiniCoeffs <- data.frame(Rule = c("A", "B", "C", "D", "E"),
                                 MeanPercentGiniChange = c(mean(SV1_RuleA_GiniCoeffs$PercentGiniChange),
                                                           mean(SV1_RuleB_GiniCoeffs$PercentGiniChange),
                                                           mean(SV1_RuleC_GiniCoeffs$PercentGiniChange),
                                                           mean(SV1_RuleD_GiniCoeffs$PercentGiniChange),
                                                           mean(SV1_RuleE_GiniCoeffs$PercentGiniChange)))

SV1_MeanGiniCoeffs <- mutate(SV1_MeanGiniCoeffs, MeanPercentGiniChange = round(MeanPercentGiniChange, digits = 3)*100)

inequalityPlot <- ggplot(data = SV1_MeanGiniCoeffs, aes(x=Rule, y = MeanPercentGiniChange)) +
  geom_bar(stat="identity") + scale_y_continuous(limits = c(-30,10), breaks = seq(-30, 10, by = 5)) +
  xlab("Bonus Vote - Rules") + ylab("Mean Percent Gini Change") +
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("SV1 - Mean Percent Gini Change")
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/SV - Inequality/SV1 - Mean Percent Gini Change.png", height = 5.5, width = 5.5)


# SV1 -  Distribution of Percent Gini Changes

SV1_GiniCoeffs_AllRules <- bind_rows(SV1_RuleA_GiniCoeffs,SV1_RuleB_GiniCoeffs,SV1_RuleC_GiniCoeffs,SV1_RuleD_GiniCoeffs,SV1_RuleE_GiniCoeffs) %>%
  mutate(PercentGiniChange = PercentGiniChange * 100)
  

inequalityPlot <- ggplot(data = SV1_GiniCoeffs_AllRules, aes(x = PercentGiniChange, colour = Rule)) +
  geom_freqpoly(binwidth = 3) + xlab("Percent Gini Change") + ylab("Absolute Count") +
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("SV1 - Smoothed Density of Percent Gini Change (per Rule)")
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/SV - Inequality/SV1 - Smoothed Density of Percent Gini Changes.png", height = 4.5, width = 8.5)






#--- SV2: Calculating Gini Coefficient, percent gini change, & freq of gini decline 

# Calculate Gini coeffs for all sims and all rules

SV2_RuleA_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "A", vote = "SV2", util_NSV = SV_IndivUtility_NSV_A, util_SV = SV_IndivUtility_SV2_A) 
SV2_RuleB_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "B", vote = "SV2", util_NSV = SV_IndivUtility_NSV_B, util_SV = SV_IndivUtility_SV2_B)
SV2_RuleC_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "C", vote = "SV2", util_NSV = SV_IndivUtility_NSV_C, util_SV = SV_IndivUtility_SV2_C)
SV2_RuleD_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "D", vote = "SV2", util_NSV = SV_IndivUtility_NSV_D, util_SV = SV_IndivUtility_SV2_D)
SV2_RuleE_GiniCoeffs <- getGiniCoeff(sampleData = SV_AllRules, rule = "E", vote = "SV2", util_NSV = SV_IndivUtility_NSV_E, util_SV = SV_IndivUtility_SV2_E)


#--- Modified ---#

# SV2_RuleA_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV2_RuleA_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# SV2_RuleB_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV2_RuleB_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# SV2_RuleC_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV2_RuleC_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# SV2_RuleD_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV2_RuleD_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# SV2_RuleE_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/Ineq - SV2_RuleE_GiniCoeffs.csv", stringsAsFactors = F, header = T)


# Gini of Average Realized Utilities (over all 10k sims)

SV2_AllRules_GiniCoeffs <- bind_rows(SV2_RuleA_GiniCoeffs, SV2_RuleB_GiniCoeffs, SV2_RuleC_GiniCoeffs, SV2_RuleD_GiniCoeffs, SV2_RuleE_GiniCoeffs)

meanGiniSummary <- bind_rows(meanGiniSummary, SV2_AllRules_GiniCoeffs %>% group_by(Rule) %>% 
                               summarize(MeanGini_Maj = mean(Gini_UtilNSV), lower95CI_Maj = getNthEmpPercentile(Gini_UtilNSV, 2.5), upper95CI_Maj = getNthEmpPercentile(Gini_UtilNSV, 97.5),
                                         MeanGini_SV = mean(Gini_UtilSV), lower95CI_SV = getNthEmpPercentile(Gini_UtilSV, 2.5), upper95CI_SV = getNthEmpPercentile(Gini_UtilSV, 97.5)) %>%
                               mutate(Sample = "SV2", AveragedOver = "All Sims"))


# Filter data so we only make plots conditional on a minority victory

SV2_RuleA_GiniCoeffs <- SV2_RuleA_GiniCoeffs %>% filter(MinVic == T) 
SV2_RuleB_GiniCoeffs <- SV2_RuleB_GiniCoeffs %>% filter(MinVic == T) 
SV2_RuleC_GiniCoeffs <- SV2_RuleC_GiniCoeffs %>% filter(MinVic == T) 
SV2_RuleD_GiniCoeffs <- SV2_RuleD_GiniCoeffs %>% filter(MinVic == T) 
SV2_RuleE_GiniCoeffs <- SV2_RuleE_GiniCoeffs %>% filter(MinVic == T) 


# Gini of Average Realized Utilities (over sims with min vic)

SV2_AllRules_GiniCoeffs <- bind_rows(SV2_RuleA_GiniCoeffs, SV2_RuleB_GiniCoeffs, SV2_RuleC_GiniCoeffs, SV2_RuleD_GiniCoeffs, SV2_RuleE_GiniCoeffs)

meanGiniSummary <- bind_rows(meanGiniSummary, SV2_AllRules_GiniCoeffs %>% group_by(Rule) %>% 
                               summarize(MeanGini_Maj = mean(Gini_UtilNSV), lower95CI_Maj = getNthEmpPercentile(Gini_UtilNSV, 2.5), upper95CI_Maj = getNthEmpPercentile(Gini_UtilNSV, 97.5),
                                         MeanGini_SV = mean(Gini_UtilSV), lower95CI_SV = getNthEmpPercentile(Gini_UtilSV, 2.5), upper95CI_SV = getNthEmpPercentile(Gini_UtilSV, 97.5)) %>%
                               mutate(Sample = "SV2", AveragedOver = "Miv Vic"))

write.csv(meanGiniSummary, "./SV & QV Recalibration on Pop Imm/SV - Inequality/Mean Gini Coeffs & CIs.csv", row.names = F)


#--- Modified (end) ---#



# Filter data so we only make plots conditional on a minority victory

SV2_RuleA_GiniCoeffs <- SV2_RuleA_GiniCoeffs %>% filter(MinVic == T) 
SV2_RuleB_GiniCoeffs <- SV2_RuleB_GiniCoeffs %>% filter(MinVic == T) 
SV2_RuleC_GiniCoeffs <- SV2_RuleC_GiniCoeffs %>% filter(MinVic == T) 
SV2_RuleD_GiniCoeffs <- SV2_RuleD_GiniCoeffs %>% filter(MinVic == T) 
SV2_RuleE_GiniCoeffs <- SV2_RuleE_GiniCoeffs %>% filter(MinVic == T) 

# SV2 - Summary of Gini Changes

SV2_AllRules_GiniCoeffs <- bind_rows(SV2_RuleA_GiniCoeffs, SV2_RuleB_GiniCoeffs, SV2_RuleC_GiniCoeffs, SV2_RuleD_GiniCoeffs, SV2_RuleE_GiniCoeffs)

rulesSummary <- SV2_AllRules_GiniCoeffs %>% group_by(Rule, ChangeDirection) %>% summarize(Count = n()) %>% mutate(Prop = round(Count/sum(Count), digits = 4))

stargazer(rulesSummary, 
          title = "Summary of Gini Changes by Rule (SV2)", type = "text", digits = 3, summary = FALSE, rownames = FALSE, 
          out = "./SV & QV Recalibration on Pop Imm/SV - Inequality/Table - Summary of Gini Changes by Rule (SV2).txt")

# SV2 - Dot Plot of Frequency of Gini Declines per Rule in (Fig 8 - top) 
FreqGiniDeclineA <- group_by(SV2_RuleA_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineB <- group_by(SV2_RuleB_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineC <- group_by(SV2_RuleC_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineD <- group_by(SV2_RuleD_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineE <- group_by(SV2_RuleE_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)

SV2_FreqGiniDeclines <- data.frame(Rule = c("A", "B", "C", "D", "E"), 
                                   Freq = c(FreqGiniDeclineA$Freq[1], FreqGiniDeclineB$Freq[1], FreqGiniDeclineC$Freq[1], 
                                            FreqGiniDeclineD$Freq[1], FreqGiniDeclineE$Freq[1]))
SV2_FreqGiniDeclines <- mutate(SV2_FreqGiniDeclines, Freq = round(Freq, digits = 3)*100)

inequalityPlot <- ggplot(data = SV2_FreqGiniDeclines, aes(x=Rule, y = Freq)) + geom_point() + 
  geom_text(aes(label=Freq), hjust = .5, vjust = 2) + ggtitle("SV2 - Percent of Gini Declines by Voting Rule") + 
  scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Freq. of Gini Declines (%)") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12))
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/SV - Inequality/SV2 - Percent of Gini Declines by Voting Rule.png", height = 5.5, width = 5.5)


# SV2 - Bar Plot of Mean Percent Gini Change per Rule in (Fig 8 - bottom) 
SV2_MeanGiniCoeffs <- data.frame(Rule = c("A", "B", "C", "D", "E"),
                                 MeanPercentGiniChange = c(mean(SV2_RuleA_GiniCoeffs$PercentGiniChange),
                                                           mean(SV2_RuleB_GiniCoeffs$PercentGiniChange),
                                                           mean(SV2_RuleC_GiniCoeffs$PercentGiniChange),
                                                           mean(SV2_RuleD_GiniCoeffs$PercentGiniChange),
                                                           mean(SV2_RuleE_GiniCoeffs$PercentGiniChange)))

SV2_MeanGiniCoeffs <- mutate(SV2_MeanGiniCoeffs, MeanPercentGiniChange = round(MeanPercentGiniChange, digits = 3)*100)

inequalityPlot <- ggplot(data = SV2_MeanGiniCoeffs, aes(x=Rule, y = MeanPercentGiniChange)) +
  geom_bar(stat="identity") + scale_y_continuous(limits = c(-30,10), breaks = seq(-30, 10, by = 5)) +
  xlab("Bonus Vote - Rules") + ylab("Mean Percent Gini Change") +
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("SV2 - Mean Percent Gini Change")
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/SV - Inequality/SV2 - Mean Percent Gini Change.png", height = 5.5, width = 5.5)


# SV2 -  Distribution of Percent Gini Changes

SV2_GiniCoeffs_AllRules <- bind_rows(SV2_RuleA_GiniCoeffs,SV2_RuleB_GiniCoeffs,SV2_RuleC_GiniCoeffs,SV2_RuleD_GiniCoeffs,SV2_RuleE_GiniCoeffs) %>%
  mutate(PercentGiniChange = PercentGiniChange * 100)


inequalityPlot <- ggplot(data = SV2_GiniCoeffs_AllRules, aes(x = PercentGiniChange, colour = Rule)) +
  geom_freqpoly(binwidth = 3) + xlab("Percent Gini Change") + ylab("Absolute Count") +
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("SV2 - Smoothed Density of Percent Gini Change (per Rule)")
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/SV - Inequality/SV2 - Smoothed Density of Percent Gini Changes.png", height = 4.5, width = 8.5)


